/*
	a)
	In this file the interface for the gradient functions are declared.
	Watershed is more effective if it is performed on a gradient of an image.
	b)
	The file gives a implementation of the immersion watershed transform as proposed by:
	- Roerdink & Meijster 
	- Luc & Soille
	with an additional function pointer to handle a passed image with markers.
*/

#ifndef _watershed_H_
  #define _watershed_H_

#include "support.h"  // templated conversion function, Image-sorting, templeated histogram
					  // this header also includes "ip.h" with functions for image processing in general, e.g. pixel-class
					  // see ip.h for more information

/* 
 what gradient should be applied?

 in:     give the input image
 otp:    chose gradient which should be applied
 param:  SE size for MM or sigma for canny-edge
 param2: sigma for gaussian blurring (scale-space)
           - the second parameter is only used if scale-space-application is desired.
*/
int applygradient(Image<short> &in, unsigned int opt=1, int param=1,int param2 = 20);

// normal MM 'grad'
void apply_grad(Image<short> &in,int s);

// eroded 'grad'
void apply_grad_e(Image<short> &in, int s);

// multi-scale grad
void apply_mg(Image<short> &in, int s);
		
// canny-edge detection
void apply_canny(Image<short> &in, int s);


	/*

		@Param in:      input image, on which watershed teransform will be performed
		@Param output:  resulting image
		@Param mark:    image containing one-pixel white inicators for every marker passed on an copy image of the input image
	
	the watershed transform segments a picture into objects.
	the watershed transform transforms an input image into output image. the output image will contain path connected regions for every local minimum.
	the function should be called on gradient images of the original source images.
	the given implementation orients on the immersion watershed analogy, proposed 
		by "Luc Vincent & Pierre Soille":
		following changes in respect to the original algorithms are implemented:
	-use of markers
	-faster(!) and more robust image sorting, via new sorting image (improves time by several seconds)
	-improvements as proposed by "Roerdink & Meijster"

	the function is overdocumented, since it is the central function.
	it should also work with real float images, which was not tested yet

	*/

template <typename T>
void watershed(const pip::Image<T> &in, pip::Image<T> &output, const pip::Image<T> mark){
	//init Images
	pip::Image<float> out(in.dimension());   // temporary image to save intermediate results
	pip::Image<T>     dist(in.dimension());  // temporary image to save distances (dist max within T-range)
	pip::Image<float> bins(in.dimension());  // temporary image to save a sorting structure, which can also be used on float and short images, without having much time issues
	out.fill(-1); // -1 is an init value
	dist.fill(0); // likewise 0 is the init value here
	//init magic numbers
	const float  WSHED(0),MASK(-2); // const values to be reused
	float curr_dist(0),label(0);    // frequently used values: label reports to which minima a pixel belongs ; distance is reporting how far away the 'next' value is. 
	bool marker_given(mark.dim(0)>1);
	//init STD	
	std::queue<Pixel> fifo;         // first in first out data structure: the basis of an breath first search, as to be performed by this algorithm
	std::vector<std::vector<Pixel> > levels; // here the pixels of the image will be sorted into a 'random access structure' so the algorithm may perform faster
	// see next function call
	// fast image sort:
	// this function sorts the image 'in' into 'levels'.
	// bins is an image, which remembers, which pixel is in which level of 'levels'
	// this is important, since there are much more levels compared to gray-values possible 

	FIsort(in,bins,levels); // the functions internal parameter 'number_of_bins' has profound impact on the performance speed!

	std::vector<Pixel> neigh;         //neigh holds the neighbours currently of interrest
	std::vector<Pixel>::iterator i,j; // iterators over neighbours
	Pixel p,p_;                       // current pixel of interest
	//marker-pixel to mark end of 'run' at level 'h'
	Pixel marker(-1,-1,-1); // marker.negative()==true
	for(unsigned int h=0; h < levels.size();++h){  // begin flooding of immersion process
		for(i= levels[h].begin(); i <  levels[h].end(); ++i){  // h points at given level
			if(i->get(bins)==(h)){                             // if pixel is at current level 
				i->set(out,MASK);							   // mask pixel: it will be reexaminated soon ...
			    i->getNeighbors(in,neigh);					   // find all neighbours of this pixel ( canditates for finding relation to nearby extrema )
				for(j = neigh.begin(); j < neigh.end(); ++j){  // find amoung all neighbours
					if(   j->get(out) >= 0){                   // all pixel, which are not watershed yet, nor labeled
						i->set(dist,T(1));					   // their distacne to the current level is one
						fifo.push(*i);						   // if we found such a pixel in the neighbourhood, push the original one into the search fifo
					}
				}
			}
		}
		curr_dist=1;											// init 
		fifo.push(marker);										// push a marker to remember this 'run' of current distance
		while(true){											// growing the current basin(s)
			p=fifo.front();										// get the front of the fifo to continue search
			fifo.pop();											// remove the element
			if(p.negative()){									// if this is end of the current 'run'
				if(fifo.empty()){								// break the current run
				   break;
				}else{											// push a newer marker for next 'run' ...
					fifo.push(marker);							
					++curr_dist;								// with higher distance
					p=fifo.front();								// and proceed with next pixel, of cause
					fifo.pop();
				}
			}
			p.getNeighbors(in,neigh);							//if not 'end of run' look at neighbours
			for(j = neigh.begin(); j < neigh.end(); ++j){		//if the neighbours are lower than current 'dist'
				if(   j->get(dist) < curr_dist
					&& j->get(out) >= 0){						// and part of basin or watershed
						if(j->get(out) > 0){
							if(   p.get(out)==MASK				// and the original pixel was at current 'h' or part of a watershed
								|| p.get(out)==WSHED){
									p.set(out,j->get(out));		// this pixel is part of it's neighbours basin 
							}else if(j->get(out)!= p.get(out)){ 
								p.set(out,WSHED);				// else it may be a watershed
							}
						}else if(p.get(out)==MASK){				
							p.set(out,WSHED);					//...
						}
				}else if(   j->get(out)==MASK					// iuf the pixel is a plateau pixel
						 && j->get(dist)==0){					
							 j->set(dist,T(curr_dist+1));		// increase the distance ( higher the more distance to border of plateu
							 fifo.push(*j);						
				}
			 }
		}
		// minima identification - with respect to marker
		for(int z=0; z < in.dim(2); ++z){						// search new minima
			for(int y=0; y < in.dim(1); ++y){					
				for(int x=0; x < in.dim(0); ++x){
					if(bins(x,y,z)==h){							// if image is in the current level (bin)
						dist(x,y,z)=0;							// and on the border of the region/plateau
						if(    out(x,y,z)==MASK					// if this is a minimum (marked? by marker)
							&& ( !marker_given || mark(x,y,z)==1)){    // a functionpointer would have many arguments... since: lets use a hack?  marker in use? no yes? where are they
							++label;							// a new labele is needed for the new minimum
							fifo.push(Pixel(x,y,z));			
							out(x,y,z)=T(label);				// this minimum is the first to be assigned to its label
							while(!fifo.empty()){				// as long as there are pixel to be searched
								p=fifo.front();
								fifo.pop();
								p.getNeighbors(in,neigh);
								for(j = neigh.begin(); j < neigh.end(); ++j){ // search neighbours
									if(j->get(out)==MASK){					  // ... which were not processed
										fifo.push(*j);						  // and put them into ...
										j->set(out,label);					  // the basin of this minimum
									}
								}
							}
						}
					}
				}
			}
		}
	}	
	convertAllToShort(out,output);
	return;
}

#endif